Localized and Cellular Patterns in a Vibrated Granular Layer 
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We propose a phenomenological model for pattern formation in a vertically vibrated layer of gran- 
ular material. This model exhibits a variety of stable cellular patterns including standing rolls and 
squares as well as localized objects (oscillons and worms), similar to recent experimental observa- 
tions(Umbanhowar et al., 1996). The model is an amplitude equation for the parametrical instability 
coupled to the mass conservation law. The structure and dynamics of the solutions resemble closely 
the properties of localized and cellular patterns observed in the experiments. 

PACS: 47.54.+r, 47.35.+i 



Granular materials exhibit a unique mixture of prop- 
erties of both liquids and solids jjj . Intensive theoretical, 
numerical and experimental studies of granular systems 
revealed a wide variety of new phenomena typical for 
granular systems, such as clustering and inelastic col- 
lapse random force chains ||, granular convection 
j| and cellular patterns in vibrated layers |§-j7|. The 
very complicated rheology of granular media makes their 
theoretical analysis extremely difficult. Unlike fluid dy- 
namics, there is no reliable continuum description of a 
granular system applicable in a wide range of conditions. 
The widely used approach is a straightforward simulation 
of many interacting particles in a gravity field [|]|]] . 

Vibrated granular systems often manifest fluid-like be- 
havior which resembles similar phenomena in conven- 
tional liquids. Recent experimental studies of vertically 
vibrated granular systems demonstrated a rich va- 
riety of collective behavior ranging from standing waves, 
hexagons, squares to localized objects (particle-like os- 
cillons |(| and one-dimensional worms, see §). Cellular 
patterns are remarkably similar to Faraday waves in flu- 
ids Jl(| which have recently been a subject of intensive 
research (see, e.g. jy]]). The difference, however, is that 
the primary bifurcation to square patterns and oscillons 
is hysteretic: patterns disappear at less magnitude of the 
plate vibrations than they first appear. The localized ob- 
jects, oscillating at half the frequency of plate vibrations 
S7/2 on the background of a flat surface oscillating at f2, 
are observed in slightly subcritical parameter region for 
cellular patterns. Although some features of collective 
behavior of vibrated granular systems were reproduced 
in simulations of a large ensemble of inelastically inter- 
acting particles |^] , the macroscopic understanding of the 
dynamics is lacking. In this Letter we propose a simple 
continuum model exhibiting phenomenology remarkably 
similar to the experimental observations. Due to its sim- 
plicity, it is amenable to a comprehensive analysis. Al- 
though this model is not derived from the corresponding 
microscopic equations of granular systems, it is instruc- 



tive for interpretation of experimental data and may yield 
testable predictions. 

The model consists of an amplitude equation for the 
order parameter ip coupled to a conservation law for an 
average mass of granular material per unit area (or a 
local averaged thickness of the layer): 

d t ip = lv b* - (1 - iuj)ip + (1 + ib)V 2 yj - \tf;\ 2 ip - pyj (1) 
d t p = aV ■ ( P V\vh\ 2 ) + 13V 2 p (2) 

Eq.(|l]) without the last term is a popular model for 
the parametric instability in oscillating liquid layer (see 
M,l3). The order parameter ijj(x,y,t) characterizes lo- 
cal complex amplitude of the particle oscillations at the 
frequency lu = Q/2. Linear terms in this equation can be 
derived from the dispersion relation for parametrically 
driven granular waves, expanded near frequency u and 
corresponding wavenumber k (here k = yJToJb, parame- 
ter b must chosen to reproduce the correct wavenumber 
at a given frequency). The term jtp* provides paramet- 
ric driving and leads to the excitation of standing waves. 
The term jV'PV' phenomenologically accounts for the non- 
linear saturation of oscillations provided in granular ma- 
terials by restitution. The last term in Eq.(|l|) accounts 
for the coupling of the order parameter to the local av- 
erage density p. As it is observed experimentally HD, 
the threshold value of the vibration amplitude 7 for the 
parametric instability depends on the mean layer thick- 
ness [g| due to an increase of internal energy dissipation in 
thicker layers. Although one should generally expect this 
term to be f(p)ip with f(p) saturating at large p (when 
the thickness of the layer is larger than the scale of typi- 
cal perturbations), we hereafter limit ourselves with the 
simplest form f(p) — p corresponding to relatively thin 
layers (proportionality constant can be omitted after ap- 
propriate scaling of p). 

Eq.(^J) describes the conservation of the granular ma- 
terial where p is a mass of granular material per unit 
square averaged over period of vibrations. Two differ- 
ent physical mechanisms contribute to the in-plane mass 
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flux. The first term in (g) reflects the average particle 
drift due to the gradient of magnitude of high-frequency 
oscillations. On average, particles try to "escape" from 
regions of large fluctuations, the effect analogous to for- 
mation of so called Chladni figures [jl2j . The second term 
describes diffusive relaxation of the inhomogeneous mass 
distribution. We expect the effective "granular tempera- 
ture" and corresponding diffusion constant /3 to be pro- 
portional to the energy of plate vibrations, then for non- 
vibrating plate the diffusion constant turns to zero and 
an arbitrary pattern "freezes" . 

Let us briefly discuss the linear stability properties of 
Eqs. The trivial state -0 = 0, p = po corresponding 

to a flat layer oscillating at the driving frequency f2 be- 
comes unstable at 7 2 = j 2 = {uj + 6(1 + p )) 2 /(l + b 2 ) for 
cub > 1 + po with respect to a periodic perturbation with 
the wavenumber k c given by k 2 — (cub — 1 — po) /(I + b 2 ). 
For cob < 1 + po spatially uniform perturbations with 
k c = become unstable first and the vibration threshold 
is 7? = 1 + p + uj 2 . Due to the rotational invariance of 
(|}g), waves with all directions grow simultaneously. 

Above the threshold, the nonlinear terms in Eqs.(|l|- 
||) saturate the exponential growth of perturbations and 
provide pattern selection. The problem of pattern selec- 
tion requires careful analysis of patterns with different 
symmetries. For p =const Eqs. reduce to a single 

equation for which it is known that rolls are the only sta- 
ble cellular pattern above onset . This has been a se- 
rious shortcoming of this model since square patterns are 
frequently observed in Faraday experiments |p|,pT|. Usu- 
ally pattern selection depends sensitively on the choice of 
nonlinearity in the model, and some tweaking with non- 
local nonlinearities could remedy this problem (see e.g. 
p5[). It turns out that within our combined model with 
p being a dynamical variable it is not needed, squares and 
rolls emerge naturally in different parameter regions. 

Consider stability of simplest cellular patterns (rolls 
and squares) in the framework of Eqs. (0-§) close to 
the threshold of parametric instability. Although the 
analysis is formally valid for arbitrary rhombic pattern, 
we restrict ourself by the squares since they are typi- 
cally preferred by the symmetry. Within the framework 
of weakly-nonlinear analysis (small supercriticality e de- 
fined below) cellular patterns are described by the fol- 
lowing solution to Eq. (fil): 



V> = (Asm(k c x) + B sm(k c y)) e 40 + w 



(3) 



where A(t), B(t) are the real amplitudes of two (orthog- 
onal) standing waves, phase (f> =const is given by so- 
lution of linearized problem, and if is a correction to 
the solution which we demand to be small at e — > 0. 
Near the threshold the density p is enslaved to |-0P, 
and follows the quasi-stationary solution of Eq. (j^) 
p = po(t) expj— r)\ip\ ], r\ — a/j3. The function po(t) 
can be found from the condition of total mass conser- 
vation S~ x J pdxdy = n = const, S is the total area. 



Substituting p into Eq. (|l]) and performing standard or- 
thogonalization procedure to keep w small, we obtain 
the following equations for A(t),B(t) (for simplicity we 
retain nonlinearity only in two first orders) 



A = A 
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Equation for B(t) is obtained by the permutation A <^=> 
B. The supercriticality parameter is given by e = 
7c(7 ~ 7c)/(l + P + k 2 c ). |f§. It follows from Eq.(|) 
that hysteretic transition to squares (A = B) occurs if 
pi] > 9/5 and stripes (A ^ 0, B = or A = 0, B ^ 0) 
exhibit subcritical bifurcation at pn > 3. 

In the supercritical case we can drop last two terms in 
Eq.(^). It is easy to verify that for e — > +0 the square 
pattern (A = B) is stable for pr] > 1 and unstable oth- 
erwise. Rolls, in the limit e — > +0, are stable for prj < 1 
and unstable otherwise. For larger e > the higher order 
terms in Eq. (^) become important, and squares are sta- 
ble at Wpen 2 < 3(10e?7 — e 2 rj 2 — 9) and rolls are stable at 
4/xer) 2 > 3(en — e 2 n 2 — 3). In the subcritical case r\p > 9/5 
squares are stable in their entire basin of existence given 
by the condition 48perj 2 + (5yp — 9) 2 = 0. The phase 
diagram for roll and squares in shown in Fig. [l]. It is 
qualitatively consistent with the experimental observa- 
tions of transition from rolls to squares with decreasing 
the driving frequency if one assumes that r\p, decreases 
with increase of frequency. One expects that the rela- 
tive effect of particle drift from regions of intense fluc- 
tuations characterized by parameter n diminishes with 
the increase of u> since characteristic vertical scale of the 
layer involved in oscillations becomes smaller at high fre- 
quencies. At large positive e there is a bistable region 
where rolls and squares co-exist, also in agreement with 
experiments ||). It should be noted however, that the 
phase diagram of Eq. (|J) exactly corresponds to the origi- 
nal model (]l|-^[) only in the limit of small amplitude, oth- 
erwise it represents only a Galerkin approximation of the 
exact solution. Still, our numerical simulations agreed 
fairly well with these stability limits. 

Now we consider the localized solutions. In experi- 
ments H oscillons appear slightly below the threshold 
of the parametric instability for cellular patterns. We 
also found stationary localized axisymmetric solutions to 
Eqs. (fiJII) in weakly subcritical region (dots in Fig. |l| 
correspond to stable localized solutions found for var- 
ious combinations of parameters). Fig. [2] shows the 
radial structure of the order parameter ip and corre- 
sponding distribution of p for such solution. This so- 
lution corresponds to a dip in the average mass dis- 
tribution p = po exp(— 77|i/>| 2 ). Due to the symmetry 
ip — > —ip, oscillons of opposite polarities may coex- 
ist in this system. It has oscillating tails at r 3> 1, 
tp(r) oc r -1 / 2 exp(pr) with the (complex) exponent p 
given by p 2 = —k 2 + (j 2 — 7 2 )/(l + b 2 ). It explains 
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that peaks and craters with elevated periphery replace 
each other on consecutive cycles of plate vibrations. 

We studied the linear stability of these localized so- 
lution with respect to axisymmetric (usually the most 
dangerous) perturbations. Some of the results are pre- 
sented in Fig. ^, where we show largest eigenvalue of 
linearized system A, 1^(0)1 and the mass deficit m — 
27r/9o(/ ° rexp[— -q\tp\ 2 ]dr — 1). The stability region is 
limited both at large and small 7 in accord with exper- 
iments. At the edges of the stable region, 7 c i,2> stable 
solution corresponding to the oscillon annihilates with 
other unstable localized solutions. 

The interaction of two oscillons can be considered in 
the spirit of Ref. Since the asymptotic behavior 

of ip for the oscillon is oscillatory, the interaction force 
F ~ Reexp(pr) is oscillatory, too. It is natural to expect 
a variety of bound states. A numerically found bound 
states of two oscillons with opposite phases is shown in 
Fig.^a, and a bound state of four oscillons (one positive 
surrounded by three negative) is shown in Fig.|]b. As in 
the experiment |(| we found that bound states with coor- 
dination numbers higher than three are unstable. There 
also exist a stable bound state of two like-phased oscil- 
lons, but the equilibrium distance between the oscillons 
is substantially larger than for oppositely-phased pairs, 
resulting in much weaker binding. Small "granular noise" 
probably unbinds such weakly coupled pairs. 

We studied the nonlinear evolution of oscillons beyond 
their region of stability in numerical simulations of Eq. 
(|I|j^). Once the lower bound 7 cl is passed, the oscillon 
rapidly decays towards a trivial state ip = 0. Increasing 
7 above 7 C 2 initial spreading lead to a range of differ- 
ent scenarios depending on other parameters r),^, u,a. 
For small r\[i (supercritical transition) oscillons produce 
a sequence of concentric rolls. Depending on r]fj, they ei- 
ther remain rolls or break to produce a disordered square 
pattern. At larger rj [i oscillon produces other oscillons 
on its periphery as seen in experiments |J. It turns 
out, surprisingly, that following oscillons do not appear 
uniformly around the center but rather organize them- 
selves in chains, resembling worm-like patterns (Fig.^Jc), 
cf. photo by P. Umbanhowar published in M. We pro- 
pose explanation of this effect by the conservation of total 
mass. Oscillons push the granular material on their pe- 
riphery. Since the excessive mass from oscillons in a chain 
is re-distributed by the diffusion, it spreads more rapidly 
near the tip of the chain. Therefore, the next oscillon 
will likely appear there, and the tip advances (compare 
with diffusion-limited growth Eq] ) . This process will con- 
tinue until the average density in the surrounding flat re- 
gions becomes so high that the creation of new oscillons is 
halted (the threshold for oscillon stability j c2 increases 
with the average density). Our simulations show that 
for values of 7 slightly above "f C 2 even after a long time, 
oscillons do not fill entire area (see Fig. |^d). 

Let us now return briefly to the problem of selection 



of cellular patterns. In the domain of oscillon stabil- 
ity these patterns can be considered as a periodic lattice 
of weakly-coupled oscillons. As it is well-known from 
the solid-state physics, the lowest energy configuration 
of like-charged objects is a hexagonal lattice (compare 
with Abrikosov lattice). In contrast, for alternatively- 
charged particles the optimal configuration is a square 
lattice where a positively-charged particle is surrounded 
by four negatively-charged (anti-ferromagnetic lattice). 
Eqs. (|l|-|^) preserve the symmetry ip — > —ip and therefore 
square patterns formed by both positive and negative 
oscillons dominates (see Fig|I|d). Once this symmetry is 
broken (as in two frequency driving), patterns are formed 
by like-phased oscillons, and the hexagonal symmetry can 
be expected. The selection of the pattern in supercriti- 
cal case where oscillons are unstable is more subtle. The 
advantage of the square pattern over rhombi cannot be 
determined in the lowest order of our weakly-nonlinear 
perturbation theory, which happens to be insensitive to 
the angle between two standing waves. Our arguments 
related to optimal packing of localized structures sug- 
gests that the unique selection of square patterns should 
occur in higher orders of e. 

We have shown on the basis of phenomenological model 
that the constraint of mass conservation plays a crucial 
role in pattern formation in vibrated granular materials. 
It leads to a subcritical bifurcation towards planar cel- 
lular patterns which naturally include squares in a wide 
range of parameters. The model also reproduces stable 
localized particle-like excitations (oscillons) and chains of 
oscillons or short rolls sometimes called worms. The pa- 
rameters of the model can be estimated experimentally 
or from molecular dynamics. We can speculate that our 
results are also relevant for fluids, where square patterns 
are ubiquitous and localized objects were observed re- 
cently fll9fl . Parametric waves on a fluid surface induce 
mean surface displacement which must obey a conserva- 
tion law. Thus, we expect that a coupled set of equations 
for the order parameter and a mean displacement may 
serve as a paradigm model for this system. Although for 
fluids the knowledge of the Navier-Stokes equations al- 
lows for the direct stability analysis of cellular patterns 
, finding a relevant order-parameter model would be 
useful for studies of more complicated pattern forma- 
tion in this rich experimental system. Worm-like struc- 
tures have also been recently observed in electroconvec- 
tion |2l| . It is plausible to assume that as in a granular 
system, this effect can also be interpreted as the growth 
in a system exhibiting a first-order transition to a cellular 
structure and controlled by a slowly diffusing field. 
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FIG. 1. Phase diagram for square and roll patterns, 
weakly-nonlinear theory. Black dots indicate stable oscillons 
seen in numerical experiments 
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FIG. 2. Radial structure of the stable oscillon for 7 = 1.8, 
H = 0.527, 6 = 2, w = a = l, 77 = 2.78 
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FIG. 3. Largest eigenvalue A (a), at r = (b) and mass 
deficit m (c) for stable (solid line) and unstable (dashed line) 
oscillons, n = 0.527, b — 2,lu — a — 1, n =5/7. 
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FIG. 4. Gray-coded images of Re-(/> (black corresponds 
to maximum, white to minimum) from simulations of 
Eqs. (Q-(^), (a) bound state of oppositely-phased oscillons, 
a = u = 1,6 = 2,77 = 2.78,/j = 0.527,7 = 1.8, size L = 40; 
(b) Triangular bound state, same parameters; (c) worm-like 
structure produced by a single oscillon in the center, a = 1, 
u = b = 2, 7 = 2.245,77 = 4.38,/* = 0.525, L = 100; (d) square 
lattice, w = a = 1, 7 = 1.84, /j = 0.52, 77 = 2.72, L = 100. 
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